#### Figure 4

# Clear
rm(list=ls())

# required packages
library(dplyr)
library(haven)


# Working directory
setwd("~/Dropbox/Paraguay/1_data/7_replication_dataverse")


### (a) Dictatorship (1954-1988) ####
# Create aggregate dataset for dictatorship period
dict <- read_dta("data/main.dta") %>%
  # Filter by period
  dplyr::filter(year < 1989) %>%
  # Paste together families and municipalities for grouping
  dplyr::mutate(family.muni = paste(lastname, district, sep = "_")) %>%
  # Select: family-municipality, eigenvector centrality, number of land grants & hectares
  dplyr::select(family.muni, eigenvector, tmh_n, tmh_ha) %>%
  # Group by family-municipality
  dplyr::group_by(family.muni) %>%
  # Aggregate: land grants & eigenvector
  dplyr::summarise(benef.fam = sum(tmh_n), total.has = sum(tmh_ha), eigenvector = mean(eigenvector)) %>%
  # Dichotomize land grants  
  dplyr::mutate(benef.fam = if_else(benef.fam >= 1, 1, 0)) %>%
  # Break eigenvector into deciles  
  mutate(decile = cut(.$eigenvector, breaks = quantile(.$eigenvector, seq(0, 1, .1), na.rm = TRUE))) %>%  
  # Now re-group by deciles
  dplyr::group_by(decile) %>%
  # Proportion of elite families that got at least a land grant
  dplyr::summarise(benef.fam = sum(benef.fam) / 16735, eigenvector = mean(eigenvector)) %>%
  # Drop NAs
  dplyr::filter(!is.na(decile) & !is.na(benef.fam) & !is.na(eigenvector)) %>%
# Enumerate deciles and calculate CIs  
  dplyr::mutate(Q = seq(1, 10, 1)) %>%
  dplyr::mutate(se_bf = sd(benef.fam) / sqrt(10)) %>%
  dplyr::mutate(upper_bf = benef.fam + (1.96*se_bf), lower_bf = benef.fam - (1.96*se_bf))

# Plot dictatorship
pdf('figures/figure4a.pdf', width = 10, height = 7)# Create base plot
plot(benef.fam ~ Q, data = dict,
     type = "l",
     lwd = 1,
     cex.lab = 1.15,
     axes = FALSE,
     ylim = c(0,.025),
     xlab = "Network Centrality (deciles)",
     ylab = "Proportion of Families",
     frame.plot = FALSE)
# Shaded area of CIs
polygon(x = c(dict$Q, rev(dict$Q)), y = c(dict$upper_bf, rev(dict$lower_bf)), col = "grey80", border = NA)
# Points
points(benef.fam ~ Q, pch = 20, col = "darkred", data = dict)
# Line
lines(benef.fam ~ Q, col = "darkred", data = dict)
# Axes
axis(side = 1, at = seq(1, 10, 1), cex.axis = 1.15)
axis(side = 2, at = seq(0, 0.025, 0.005), cex.axis = 1.15)
dev.off()

### (b) Democracy (1989-2003) ####
# Create aggregate dataset for democracy period
demo <- read_dta("data/main.dta") %>%
  # Filter by period
  dplyr::filter(year > 1989) %>%
  # Paste together families and municipalities for grouping
  dplyr::mutate(family.muni = paste(lastname, district, sep = "_")) %>%
  # Select: family-municipality, eigenvector centrality, number of land grants & hectares
  dplyr::select(family.muni, eigenvector, tmh_n, tmh_ha) %>%
  # Group by family-municipality
  dplyr::group_by(family.muni) %>%
  # Aggregate: land grants & eigenvector
  dplyr::summarise(benef.fam = sum(tmh_n), total.has = sum(tmh_ha), eigenvector = mean(eigenvector)) %>%
  # Dichotomize land grants  
  dplyr::mutate(benef.fam = if_else(benef.fam > 1, 1, 0)) %>%
  # Break eigenvector into deciles  
  mutate(decile = cut(.$eigenvector, breaks = quantile(.$eigenvector, seq(0, 1, .1), na.rm = TRUE))) %>%  
  # Now re-group by deciles
  dplyr::group_by(decile) %>%
  # Proportion of elite families that got at least a land grant
  dplyr::summarise(benef.fam = sum(benef.fam) / 16735, eigenvector = mean(eigenvector)) %>%
  # Drop NAs
  dplyr::filter(!is.na(decile) & !is.na(benef.fam) & !is.na(eigenvector)) %>%
# Enumerate deciles and calculate CIs
  dplyr::mutate(Q = seq(1, 10, 1)) %>%
  dplyr::mutate(se_bf = sd(benef.fam) / sqrt(10)) %>%
  dplyr::mutate(upper_bf = benef.fam + (1.96*se_bf), lower_bf = benef.fam - (1.96*se_bf))

# Plot democracy
pdf('figures/figure4b.pdf', width = 10, height = 7)# Create base plot
plot(benef.fam ~ Q, data = demo,
     type = "l",
     lwd = 1,
     cex.lab = 1.15,
     axes = FALSE,
     ylim = c(0,.025),
     xlab = "Network Centrality (deciles)",
     ylab = "Proportion of Families",
     frame.plot = FALSE)
# Shaded area of CIs
polygon(x = c(demo$Q, rev(demo$Q)), y = c(demo$upper_bf, rev(demo$lower_bf)), col = "grey80", border = NA)
# Points
points(benef.fam ~ Q, pch = 20, col = "darkred", data = demo)
# Line
lines(benef.fam ~ Q, col = "darkred", data = demo)
# Axes
axis(side = 1, at = seq(1, 10, 1), cex.axis = 1.15)
axis(side = 2, at = seq(0, 0.025, 0.005), cex.axis = 1.15)
dev.off()




